Condensate fraction of cold gases in non-uniform external potential 
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Exact calculation of the condensate fraction in multi-dimensional inhomogeneous interacting Bose 
systems which do not possess continuous symmetries is a difficult computational problem. We have 
developed an iterative procedure which allows to calculate the condensate fraction as well as the 
corresponding eigenfunction of the one-body density matrix. We successfully validate this procedure 
in diffusion Monte Carlo simulations of a Bose gas in an optical lattice at zero temperature. We 
also discuss relation between different criteria used for testing coherence in cold Bose systems, such 
as fraction of particles that are superfluid, condensed or are in the zero-momentum state. 
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Bosc-Einstcin condensation (BEG) is a fascinating 
phenomenon in which the effects of quantum coherence 
become apparent on a macroscopic scale [Jill • Although 
the phenomenon is knovifn for a long time and the ba- 
sic theoretical concepts are well established, there is no 
general method which would allow to calculate exactly 
such fundamental quantities as the condensate fraction 
Nq/N and the wave function of the condensate 0o(r), 
which can be probed in modern experiments with ultra- 
cold atoms (see, e.g., references in [l|-|3[). In order to cal- 
culate those, one should solve the eigenvalue problem for 
one-body density matrix (OBDM) pi(r;r') [ll 0] which 
is generally a very difficult task. Simple solutions can 
be obtained for weakly interacting gas at zero tempera- 
ture. In this case the condensation is almost complete, 
Nq/N « 1, and 4>o{^) can be obtained by means of the 
mean-field theory as a solution of the Gross-Pitaevskii 
equation (GPE) ^'^■^'^(r) Q. Quantum corrections to 
the mean-field predictions at zero temperature can be 
calculated within the framework of the Bogoliubov the- 
ory Q. 

Monte Garlo (MG) methods make no approximation 
to the model as compared to the perturbative techniques 
and can address the case of strong interactions. In Ref. Q 
properties of a hard-core Bose gas in a harmonic trap 
were studied by MG methods. Taking into account spher- 
ical symmetry of the confinement and making some ad- 
ditional assumptions, which are valid only for weakly in- 
teracting gas, the original problem in three dimensions 
was reduced to an effective onc-dimcnsional problem and 
then standard methods of matrix diagonalization were 
applied to calculate (/)o(r)- However, in a large number of 
experiments with ultracold atoms the external potentials 
do not possess spherical symmetry. This is the situation 
of a gas in the optical lattice 0, Q , gas in the presence 
of a disorder potential Q, etc. In such cases there is a 
non-trivial dependence of OBDM on all arguments and 
it is not possible to reduce the diagonalization problem 



to a simple matrix formulation. In this Letter we pro- 
pose a general method for OBDM with arbitrary coordi- 
nate dependencies in higher spatial dimensions. We will 
also compare the condensate fraction with other quanti- 
ties used for quantitative description of coherence such 
as the fraction of condensed particles in the momentum 
space and the superfluid fraction. 

The eigenvalue problem for the OBDM reads as 



/ Pi(r;r') 



,(r')dr'-7V,0,(r) , ^ = 0,l,, 



(1) 



with eigennumbers Ni labeled in descending order and 
eigenfunctions satisfying the orthonormality condition 



^(r)<^,(r)rfr-<5. 



In order to work out the largest eigenvalue Nq and the 
corresponding wave function (t>o(r), we use an idea stem- 
ming from the matrix analysis. When a matrix acts on 
a vector it produces a new vector which can be obtained 
from the initial one by multiplication of its components 
along the directions of the eigenvectors by the corre- 
sponding eigenvalues. If the resulting vector is renor- 
malized such that it has the same norm as the original 
one, the result of applying the matrix is a rotation of 
the vector in the direction of the eigenvector with the 
largest eigenvalue. Iterating such a rotation many times 
will eventually align the original vector with the eigen- 
vector with the largest eigenvalue. The convergence of 
the iterative procedure is very fast if one of the eigenval- 
ues is much larger than the others which is the case of a 
Bose-condensed system. 

Applying this idea to Eq. ([T]) , we come to the iterative 
procedure, where the (i + l)-th approximation for the 
wave function is determined by 

y'pi(r;r')0«(r')dr'=iV«</>r'^(r). (2) 
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The i-th approximation for the number of condensed par- 
ticles is given by 



N, 



p^iv;v')4^*ir)4\r')drdr' , (3) 



which follows from Ecj^. ([TJ. Repeating this procedure 
permits to obtain, in principle, the condensate wave func- 
tion exactly. A reasonable choice for the initial approxi- 
mation is (/ig^'' (r) = (fp^^ (r) . 

During the iterations the value of Nq is approached 
from below. This can be seen by first expand- 
ing 4>^q\v) in terms of eigenfunctions of the OBDM 
(r) = X]^o "^J '^J with the normalization condition 
Sj^o ^ ^ ^'^"^ then inserting the resulting expres- 
sion to Eq. ([3]). This leads to inequality 



N, 



OO 



i=0 



where A^i is the upper bound for the eigenvalues with 
1 = 1,2,... This proves the statement. 

The iterative procedure can be implemented in MC 
calculations. The OBDM in first quantized form is ex- 
pressed in terms of the many-body wave function of the 
ground state ^/'(R) as [HQ 

Pl(r;r')=Ar j r(R)lri=r^(R)lri=r' dV2...dVM, 

where R = {yi, . . . ,ym) is a shortcut for a point in 
3A'^-dimensional phase space. This allows to rewrite 
Eqs. (I2|), ([3]) in a form which can be interpreted in terms 
of a MC algorithm 



AT, 



N 



i/)*(R) 



'0*(r,r2,...,rAr) 



|^/'(R)rdR, 



dr 



\^(R)rdR. 



MC calculation Q produces a set of configurations 
Ri,R2, ... in the phase space, distributed according to 
the best approximation of |?/'(R)|^. The averaging of 
each of the quantities in square brackets is then equiva- 
lent to the averaging over the produced set of configura- 
tions Ri, R2, .... In the case of variational MC method. 
Metropolis sampling of the trial function -0t(R) pro- 
duces configurations distributed according to |'0t(R)| ■ 
The diffusion Monte Carlo (DMC) method samples the 
"mixed" distribution -0t(R)V'(R-)- "Pure" averages of 
A^o and 0o over the ground state wave function liAlR-) 
are approximated by extrapolation from variational and 
mixed estimators. The integral over r can be evaluated 
in a stochastic way by sampling a random point in the 
simulation box and accumulating the values of the aver- 
aged quantity. Since the OBDM is calculated as a mixed 



estimator, the result that A^q*^ is a lower bound to A^o 
[Eq. dU] should not necessarily hold. 

Sufficient condition for the existence of BEC is the off- 
diagonal long-range order (ODLRO) of the OBDM 0, i, 
IToj . If the asymptotic value 



A^, 



k=0 



lim pi(r;r') 



r— r' — >oo 



(5) 



does not vanish, there is a finite fraction of particles with 
zero momentum (k = 0). A"^^" can be used as a mea- 
sure of BEC in homogeneous and slightly inhomogeneous 
systems H, llj. However, in general the wave function 
of the state with k = is not the eigenfunction of the 
OBDM and as it follows from Eq. (g]) N^=° < No- 

Another quantity used to describe coherence in inter- 
acting quantum systems is the superfluid fraction. The 
number of atoms in the superfluid can be obtained as [l3| 
Ns = \imv^o2AF/{mv^), where AF is the increase of 
the free energy in the reference frame moving with the 
velocity v. It is interesting to compare Ng with A^o £md 

In a weak external potential V(r) with the period L in 
d spatial dimensions, A"^ 
perturbative framework as [1 



Q can be approximated within 



N 
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(6) 



where gd is an effective interaction parameter in d di- 
mensions, Ud = N/L'^ is the number of atoms per unit 
volume. For a dilute gas with spherically symmetric in- 
teraction potential, gd is proportional to the s-wave scat- 
tering length tts- V(n) is the Fourier transform of the 
external potential V{r), i.e.. 



V{n) 
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L/2 



l^d/2 



dr\ 



dvd e" 



-L/2 



-L/2 



Equation (|6| is obtained as a perturbative solution of 
the CPE [l5[ and does not take into account quantum 
fluctuations (Lee-Huang- Yang correction Analo- 
gous calculations for the superfluid fraction lead to an 
expression similar to Eq. ([5]) but with the coefficient 4/(i 
in front of the sum. The same result follows from the 
Bogoliubov theory. We note that known results for sys- 
tems with 5-correlated disorder can be reproduced 
by Eq. (jG)) after statistical averaging. 

In order to make a direct comparison between A'o, 
Nq^^, and Ns, we do numerical simulations of a Bose 
gas in an optical lattice described by the following many- 
body Hamiltonian 



H 



N 

E 

i=l 



2m 



,1), (7) 
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where Vpp{r) is a particle-particle interaction potential. 
In the MC calculations, we use the hard-sphere potential 
of the radius Og. The GPE is solved for the i5-potential 
characterized by the scattering length equal to a^. The 
external potential has the form 



y(r) = Fo^cos2 (2^^ 

q:=1 ^ ^ 



(8) 



where Al/2 is the lattice period. For simplicity we con- 
sider a quasi-two dimensional {d = 2) geometry when the 
system is subjected to such a tight harmonic oscillator 
trapping in the third dimension V{z) = nnJ^^z^ /2 that 
the energy of particles is small compared to the energy of 
the harmonic confinement hui]xo- This corresponds to re- 
cent experiments in anisotropic optical lattices, where the 
confinement in the third dimension was produced by a 
periodic potential of large amplitude [111. For this setup, 
Eq. ([ni reduces to 



N 



1-^ — 



{2Er + gdrid) 



(9) 



where En — 2h^Ti'^ / {m\\) is the recoil energy. 

In our calculations, the system parameters are cho- 
sen to remain in the superfluid part of the phase dia- 
gram. Mean-field theory in the tight-binding approxima- 
tion gives the following critical value for the superfluid- 
Mott-insulator transition 0]: 



2dJ/Ud = 2?I + 1 - 2^n{n+l) , 



(10) 



where n is the number of atoms per lattice site, which 
must be integer, J is the tunneling rate, and Ud is the 
interaction parameter. For d — 2, we get an estimate 



Ud as V TTflho \huJhoJ ^ \ TTflho V ^ho 



exp -- 

\ 



2Fn 



(11) 



where Oho = x/^i/mwho is the harmonic oscillator length. 
In order to remain in the superfluid regime, the values 
of 2dJ/Ud should be larger than that given by Eq. ([T0| . 
which leads to restrictions on the values of Vq and . 

The one-body part of the variational wave function 
iI)t{y) used in MC calculation is obtained by solving GPE 
for a single lattice period. A plausible approximation 
for the condensate orbital is 0o(r) = (fP ^ ^ {x , y)^\io{z) , 
where V'ho(-z) is the ground-state wave function of the 
harmonic oscillator. We fix the height of the optical lat- 
tice to Vb = O.Sfiwho and restrict ourselves to the case 
of one atom per lattice period {n = 1). The calculations 
are carried out for Ai/aho = 15, which corresponds to 
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the experimental setup in Ref. 

We start the numerical investigation with the case 
of a weak interaction ag/oho = 0.1. According to 
Eqs. ((TU)) . (|lll) the critical value of the lattice strength 




Figure 1. One-body density matrix pi{r) normalized to mean 
density 712 in a system of = 100 particles as obtained by 
extrapolation procedure in DMC calculation. The parame- 
ters are Vb = O.S&Jho, fls/aho = 0.1 (squares, upper curve), 
1 (circles, lower curve). 



is = 0.89 ftwho- It turns out from MC calculations 
that the condensate fraction is very large Nq/N w 0.99 
in a wide range of the strengths of the optical lattice. 
This means that (i) the system is fully condensed, (ii) 
the guess 0o(r) — (f)'^^^ {x,y)'ilJho{z) for the condensate 
orbital is indeed extremely good. Figure [1] shows the av- 
eraged OBDM pf (r) = J^J dzpi(r' + r,z;r',z), 
where r and r' are two-dimensional vectors. Its long- 
range asymptotic value gives the fraction of particles 
with zero momentum in (x — y) plane. As it is seen 
from the figure, Nq^'^/N w 0.85 which is considerably 
smaller than the condensate fraction Nq/N. We consider 
this as a convincing example that the number of par- 
ticles in the condensate should not be confounded with 
the number of particles with zero momentum. The GPE 
approach works very well in the dilute regime and pre- 
dicts the same fraction of particles with zero momentum, 
while the result of the Bogoliubov theory is slightly lower 
N^=°/N w 0.83. 

The superfluid fraction in the considered case calcu- 
lated by DMC method Ns/N = 0.75 coincides with the 
value obtained from the GPE. Perturbative Bogoliubov 
theory gives Ns/N = 0.67, refer to Eq. (P. It is inter- 
esting to note that the superfluid fraction is gradually 
reduced with increasing Vb, while the condensate frac- 
tion remains very large. A possible interpretation is that 
for such small values of the gas parameter the external 
potential effectively changes very smoothly (i.e. classi- 
cally), so the system remains well described by the GPE, 
which corresponds to having almost all particles in the 
condensate. At the same time the superfluid flow of par- 
ticles becomes blocked by the strong external field. In 
this situation, Ns can be much smaller than A^o- 

Next we study a situation when the condensate frac- 
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Figure 2. Diagonal terms x = y oi the projected condensate 
wave function (\f^(x,y) = J <j}o{x,y,z) dz obtained as an iter- 
ation ((2]) of the solution of GPE in {x, y) plane multiplied by 
the Gaussian in z direction. 



gives a fraction which is smaller than the correct conden- 
sate fraction. We propose to approximate the conden- 
sate orbital by a solution of GPE and prove that corre- 
sponding occupation number is a lower bound to Nq. We 
cheek in DMC calculations that such an approximation 
to the condensate orbital can be successfully used even 
in strongly interacting Bose gases. Numerical results are 
compared to predictions of perturbative Bogoliubov the- 
ory. 
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tion is small. To do so, we consider a strongly inter- 
acting case with Os = aho, for which Eqs. (jlOp . (fTTj) 
predict transition at = 0.32 huy^o. For this in- 
teraction the quantum fluctuations deplete the conden- 
sate in a homogeneous system by approximately 20%. 
Presence of a lattice reduces the condensate fraction to 
No/N = 0.7. The zero-momentum fraction N^=°/N 
is further diminished to about 0.45 (see Fig. [T]). GPE 
as well as Eq. (O give 0.92. Similarly, from the GPE 
we get Ns/N = 0.87, while the Bogoliubov theory pre- 
dicts a close value N^/N = 0.84. Instead, DMC result 
Ns/N = 0.6 is significantly smaller. This large discrep- 
ancy between MC and GPE is not only due to the strong 
infiuence of quantum fiuctuations but also due to the dif- 
ferent forms of the atomic interaction potential. 

Finally, we test whether the solution of the GPE repro- 
duces well the condensate orbital 4>o{r). The results are 
presented in Fig. [2] for the case of large s-wave scattering 
length cis = aho- We find that even in such a strongly 
interacting system the condensate orbital constructed as 
a product of (fP^^{v) by the Gaussian indeed turns out 
to be almost an eigenstate. On the variational level the 
Jastrow terms are responsible for suppression of the con- 
densate fraction, although in the studied case the shape 
of the VMG orbital remains almost unaffected and is very 
close to the solution of GPE. The DMC algorithm cor- 
rects the orbital and makes it less localized compared to 
the GPE prediction (see Fig. [2|) . Similar effect has been 
observed in harmonic traps where the condensate moves 
to the edges 0. 

To conclude, we have developed a procedure of obtain- 
ing the number of condensed particles Nq which is ap- 
plicable to inhomogeneous systems in higher dimensions. 
For the experimentally relevant case of the Bose gas in 
the optical lattice, we show that the frequently used cri- 
teria of number of particles with zero momentum N^^^ 
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